c   program name: MESamp.for.

c    Tests user supplied pairwise correlation structure

c     for internal consistency, performs Cholesky decomposition &

c     generates correlated sample numbers based upon rank correlation.

c     The resultant sample numbers will be used to aggregate ANWR play

c     results.

c

c    Written 4/2/98, Schuenemeyer

c

c    Input file:

c      Unit 10 - DepA.dat, Contains play number, number of plays, play 

c           names and correlation matrix

c    Output file:

c      Unit 11 - Rand.dat, contains sample index numbers

c

c    Subroutines required: jacobi.for dcholdc.for buble.for

c      & bublei.for (included)

c

      parameter(epsilon=0.001,lx=50)

      character pln*15,title*70

	real*8 wk2(50),dtri(50,50)

	real mult

      dimension rho(50,50),tri(50,50),nply(50)

     1 ,asvd(50,50),ev(50),v(50,50),wk1(lx),id(lx)

     2 ,iv(50),pln(50),u(50),au(50),x(10000,11)

	3 ,is(10000),st(10000),it(10000)

      write(*,2)

    2 format(' Program MESamp.for generates correlated sample indices')

c   open correlation matrix file

      open(10,file='DepA.dat',status='old')

	open(11,file='Rand.dat')

      do i=1,np

       id(i)=i

      end do

      nerr=0

c

c   read correlation file

	read(10,'(a70)') title

c      np is number of plays

	read(10,*) np		

	do i = 1,np

c      iv is play num, nply is num of replications, pln is play name

	 read(10,'(i2,1x,i5,1x,a15)')iv(i),nply(i),pln(i)

	end do

c   read "correlation matrix"

	do i=1,np

       read(10,*) (rho(i,j),j=1,i)

	end do

	do i=1,np

       do j=1,i 

        rho(j,i)=rho(i,j)

        asvd(i,j)=rho(i,j)

        asvd(j,i)=rho(i,j)

        tri(i,j)=0.0

        tri(j,i)=0.0

	 end do

      end do

c   compute eigenvalues; see if `corr matrix' pos semi-def

      do i=1,np

       rho(i,i)=1.0

       asvd(i,i)=1.0

       tri(i,i)=0.0

      end do

      call jacobi(asvd,np,lx,ev,v,nrot)

      do i=1,np

       wk1(i)=ev(i)

      end do

      call buble(wk1,id,np)

      write(*,15)

   15 format(/' Eigenvalues'/)

      do i=1,np

       write(*,17)i,wk1(i)

   17  format(i5,f7.3,5f7.3)

      end do

c   allows user to view eigenvalues

      pause 0

      if (wk1(1).ge.0.0) then

       write(*,19) wk1(1)

   19 format(/' Minimum eigenvalue=',f10.3/

     1 3x,'Matrix is correlation matrix-bias factor not needed')

	 do i=1,np

	  dtri(i,i)=1.0

	  do j=1,i

	   dtri(i,j)=rho(i,j)

	   dtri(j,i)=dtri(i,j)

	  end do

	 end do

	 goto 40

      end if

      bias=abs(wk1(1))+epsilon

      write(*,22)wk1(1),bias

   22 format(/' Min. eigenvalue=',f10.3,'. Bias=',f10.3)

      do i=1,np

       ev(i)=ev(i)+bias

      end do

c  Generate new adjusted correlation matrix

      do i=1,np

      do j=1,np

        do k=1,np

         tri(i,j)=tri(i,j)+v(i,k)*ev(k)*v(j,k)

        end do

      end do

      end do

c  Normalize new correlation matrix.

      xnc=tri(1,1)

      do i=2,np

       do j=1,i-1

        tri(i,j)=tri(i,j)/xnc

        tri(j,i)=tri(i,j)

	 end do

      end do

	do i=1,np

	tri(i,i)=1.0

	dtri(i,i)=1.0

	 do j=1,i

	  dtri(i,j)=tri(i,j)

	  dtri(j,i)=dtri(i,j)

	 end do

	end do

c

c    Compute Cholesky decomposition to get triangular matrix

	call dcholdc(dtri,np,lx,wk2)

      write(*,18)

   18 format(/' New Correlation Matrix-upper right - '

     1 'Cholosky-lower left '/)

	do i=1,np

	write(*,33)iv(i),(dtri(i,j),j=1,np)

   33 format(i3,10f6.3)

      end do

c   Perform check on Cholesky decomp

   40 tdmax=0.

	do i=1,np

	 do j=i,np

	   t=0

	   do k = 1,i

	    t=t+dtri(j,k)*dtri(i,k)

	   end do

	   if(i.eq.j) then

	   td=abs(t-1.)

	   else

	   td=abs(t-dtri(i,j))

	   end if

	   if(td.gt.tdmax) dtmax=td

	 end do

	end do

	if (abs(td).gt.0.001) then

	  write(*,*) ' Warning,Cholesky decomp suspect'

	  write(*,*) '   Deviations too large', td

	end if

c	   

c   Generate random correlated sample numbers

c

 	iseed=871

	call seed(iseed)

c   num is number of samples

	num=10000

	do m=1,num

	 do j=1,np

c   generates uniforms in range (-1,+1)

	  u(j)=rand()*2. - 1.

       end do

c   multiply lower triangular * u to get dependent structure

       do i=1,np

	  au(i)=0.0

	  do j=1,i

	   au(i)=au(i)+dtri(i,j)*u(j)

	  end do

	  x(m,i)=au(i)

       end do

      end do

c   rank the sample

	do i=1,np

	 do m=1,num

	  is(m)=m

	  it(m)=m

	  st(m)=x(m,i)

	 end do

	 call buble(st,is,num)

	 call bublei(is,it,num)

	 mult=1.0

	 if(nply(i).gt.10000) mult=real(nply(i))/10000.

	 do m=1,num

	  x(m,i)=int(real(it(m))*mult)

	 end do

	end do

c this is for Nig(many), play 11

      do m=1,num

	 x(m,11)=int(x(m,10)*1.5432)

	end do

c write sample index numbers

	do m=1,num

	 write(11,58)m,(x(m,j),j=1,np+1)

   58  format(i6,11f8.0)

      end do

      stop

      end

	SUBROUTINE jacobi(a,n,np,d,v,nrot)

C     Numerical Recipes Software, 1986-92

c      d is eigenvalaues of a in first n elements.

c      v contains the normalized eignevectors.

c      note d=v trans * a * v

      INTEGER n,np,nrot,NMAX

      REAL a(np,np),d(np),v(np,np)

      PARAMETER (NMAX=500)

      INTEGER i,ip,iq,j

      REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)

      do 12 ip=1,n

        do 11 iq=1,n

          v(ip,iq)=0.

11      continue

        v(ip,ip)=1.

12    continue

      do 13 ip=1,n

        b(ip)=a(ip,ip)

        d(ip)=b(ip)

        z(ip)=0.

13    continue

      nrot=0

      do 24 i=1,50

        sm=0.

        do 15 ip=1,n-1

          do 14 iq=ip+1,n

            sm=sm+abs(a(ip,iq))

14        continue

15      continue

        if(sm.eq.0.)return

        if(i.lt.4)then

          tresh=0.2*sm/n**2

        else

          tresh=0.

        endif

        do 22 ip=1,n-1

          do 21 iq=ip+1,n

            g=100.*abs(a(ip,iq))

            if((i.gt.4).and.(abs(d(ip))+

     1       g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then

              a(ip,iq)=0.

            else if(abs(a(ip,iq)).gt.tresh)then

              h=d(iq)-d(ip)

              if(abs(h)+g.eq.abs(h))then

                t=a(ip,iq)/h

              else

                theta=0.5*h/a(ip,iq)

                t=1./(abs(theta)+sqrt(1.+theta**2))

                if(theta.lt.0.)t=-t

              endif

              c=1./sqrt(1+t**2)

              s=t*c

              tau=s/(1.+c)

              h=t*a(ip,iq)

              z(ip)=z(ip)-h

              z(iq)=z(iq)+h

              d(ip)=d(ip)-h

              d(iq)=d(iq)+h

              a(ip,iq)=0.

              do 16 j=1,ip-1

                g=a(j,ip)

                h=a(j,iq)

                a(j,ip)=g-s*(h+g*tau)

                a(j,iq)=h+s*(g-h*tau)

16            continue

              do 17 j=ip+1,iq-1

                g=a(ip,j)

                h=a(j,iq)

                a(ip,j)=g-s*(h+g*tau)

                a(j,iq)=h+s*(g-h*tau)

17            continue

              do 18 j=iq+1,n

                g=a(ip,j)

                h=a(iq,j)

                a(ip,j)=g-s*(h+g*tau)

                a(iq,j)=h+s*(g-h*tau)

18            continue

              do 19 j=1,n

                g=v(j,ip)

                h=v(j,iq)

                v(j,ip)=g-s*(h+g*tau)

                v(j,iq)=h+s*(g-h*tau)

19            continue

              nrot=nrot+1

            endif

21        continue

22      continue

        do 23 ip=1,n

          b(ip)=b(ip)+z(ip)

          d(ip)=b(ip)

          z(ip)=0.

23      continue

24    continue

      pause 'too many iterations in jacobi'

      return

      END

	SUBROUTINE dcholdc(a,n,np,p)

C  Numerical Recipes Software 

c   Corrected by schuenemeyer 3/4/94.

c   Cholesky square root decomposition of n x n matrix.

c   Stored in lower triangular of a.

c   Double precision version

      INTEGER n,np

      REAL*8 a(np,np),p(n) ,sum

         

      do 13 i=1,n

        do 12 j=i,n

          sum=a(i,j)

	    if(i.gt.1) then

          do k=i-1,1,-1

            sum=sum-a(i,k)*a(j,k)

          end do

	    end if

          if(i.eq.j)then

            if(sum.le.0.) then

		   write(*,*) 'choldc failed',i,sum

	       pause 9

	      end if

            p(i)=dsqrt(sum)

            a(i,i)=p(i)

          else

            a(j,i)=sum/p(i)

          endif

12      continue

13    continue

      return

      END

	SUBROUTINE BUBLE(X,ID,N)

c   Sorts X in ascending order, carries along ID

c   X is real, N is number of observations in X.

      DIMENSION X(1),ID(1)

      KS=N

   15 KW=0

      DO 30 I=2,KS

      IF(X(I).GE.X(I-1)) GOTO 30

      TMP=X(I)

      X(I)=X(I-1)

      X(I-1)=TMP

      NTI=ID(I)

      ID(I)=ID(I-1)

      ID(I-1)=NTI

      KW=1

   30 CONTINUE

      IF(KW.EQ.0) RETURN

      KS=KS-1

      IF(KS.EQ.1) RETURN

      GOTO 15

      END

	SUBROUTINE BUBLEI(X,ID,N)

c   Sorts X in ascending order, carries along ID

c   X is integer, N is number of observations in X.

      integer x,tmp

      DIMENSION X(1),ID(1)

      KS=N

   15 KW=0

      DO 30 I=2,KS

      IF(X(I).GE.X(I-1)) GOTO 30

      TMP=X(I)

      X(I)=X(I-1)

      X(I-1)=TMP

      NTI=ID(I)

      ID(I)=ID(I-1)

      ID(I-1)=NTI

      KW=1

   30 CONTINUE

      IF(KW.EQ.0) RETURN

      KS=KS-1

      IF(KS.EQ.1) RETURN

      GOTO 15

      END

